!##############################################################!
!                                                              !
!	    2D Incompressible Navier-Stokes Equations Solver       !
!                                                              !
! NS Eqs : dQdT + dEdX +dFdY = 0                               !
! Author : He Chen                                             !
! Time   : 2018/01/25                                          !
! Method : Artificial Compressibility                          !
! Refs   : [1] Chorin, 1967                                    !
!          [2] Hoffmann, CFD                                   !
!##############################################################!


!**************************************************************!
!                                                              !
!                    VARIABLES DEFINITION                      !
!                                                              !
! XI   --- X-coordinate of mesh points in physical domain      !
! ET   --- Y-coordinate of mesh points in physical domain      !
! X    --- X-coordinate of mesh points in computational domain !
! Y    --- Y-coordinate of mesh points in computational domain !
! U    --- Velocity in X direction                             !
! V    --- Velocity in Y direction                             !
! P    --- Pressure                                            !
! dT   --- Time step                                           !
! dX   --- X increment                                         !
! dY   --- Y increment                                         !
! RE   --- Reynolds number                                     !
! EPS  --- Convergence criteria                                !
! dUdX --- Partial dirivative                                  !
! dUdY --- Partial dirivative                                  !
! dVdX --- Partial dirivative                                  !
! dVdY --- Partial dirivative                                  !
! dPdX --- Partial dirivative                                  !
! dPdY --- Partial dirivative                                  !
!**************************************************************!

PROGRAM INSES
IMPLICIT NONE

INTEGER :: I, J, M, N, ITER, MAXITER
REAL *8 :: XIMAX, ETMAX, DELT, OMEGA
REAL *8 :: RE, dX, dY, dT, EPS, TEMP, RES
REAL *8 :: dXdET, dYdXI, ALPHA, BETA, PTENSOR

REAL *8, ALLOCATABLE :: XI(:), ET(:)
REAL *8, ALLOCATABLE :: dXdXI(:), dYdET(:)
REAL *8, ALLOCATABLE :: U1(:, :), U2(:, :), U3(:, :)
REAL *8, ALLOCATABLE :: V1(:, :), V2(:, :), V3(:, :)
REAL *8, ALLOCATABLE :: P1(:, :), P2(:, :), P3(:, :)

!===========================================!
! 			PARAMETER SPECIFICATION         !
!===========================================!
RE    = 1.0D0
XIMAX = 1.0D0
ETMAX = 1.0D0
DELT  = 0.00032D0
EPS   = 1E-10

MAXITER = 1000

M  = 19
N  = 19

dX = 1.0D0 / (M - 1)
dY = 1.0D0 / (N - 1)
dT = 0.60D0 * dX * DELT ** 0.50D0

ALLOCATE(XI(M))
ALLOCATE(ET(N))

ALLOCATE(dXdXI(M))
ALLOCATE(dYdET(N))

ALLOCATE(U1(M, N))
ALLOCATE(U2(M, N))
ALLOCATE(U3(M, N))

ALLOCATE(V1(M, N))
ALLOCATE(V2(M, N))
ALLOCATE(V3(M, N))

ALLOCATE(P1(M, N))
ALLOCATE(P2(M, N))
ALLOCATE(P3(M, N))

!===========================================!
! 			   MESH GENERATION              !
!===========================================!
    !ALPHA = 0.50D0
    !BETA  = 1.050D0
    !
    !DO I = 0, M - 1
    !		   TEMP  = I * dX
    !		   
    !	       TEMP  = ((BETA + 1) / (BETA - 1)) ** ((TEMP - ALPHA) / (1 - ALPHA))
    !		   
    !	   XI(I + 1) = XIMAX * ((2 * ALPHA + BETA) * TEMP + 2 * ALPHA - BETA) /  &
    !				   ((2 * ALPHA + 1) * (1 + TEMP))
    !				   
    !	dXdXI(I + 1) = 2 * BETA * (2 * ALPHA + 1) * (1 - ALPHA) / (XIMAX * (BETA &
    !				   ** 2 - ((2 * ALPHA + 1) * XI(I + 1) / XIMAX - 2 * ALPHA)  &
    !				   ** 2) * LOG((BETA + 1) / (BETA - 1)))
    !END DO
    !
    !
    !ALPHA = 0.50D0
    !BETA  = 1.050D0
    !
    !DO J = 0, N - 1
    !		   TEMP  = J * dY
    !		   
    !	       TEMP  = ((BETA + 1) / (BETA - 1)) ** ((TEMP - ALPHA) / (1 - ALPHA))
    !		   
    !	   ET(J + 1) = ETMAX * ((2 * ALPHA + BETA) * TEMP + 2 * ALPHA - BETA) /  &
    !				   ((2 * ALPHA + 1) * (1 + TEMP))
    !				   
    !	dYdET(J + 1) = 2 * BETA * (2 * ALPHA + 1) * (1 - ALPHA) / (XIMAX * (BETA &
    !				   ** 2 - ((2 * ALPHA + 1) * ET(J + 1) / ETMAX - 2 * ALPHA)  &
    !				   ** 2) * LOG((BETA + 1) / (BETA - 1)))
    !END DO

DO I = 1, M
    XI(I) = (I - 1) * dX
END DO

DO J = 1, N
    ET(J) = (J - 1) * dY
END DO

dXdXI = 1.0D0
dYdET = 1.0D0

!===========================================!
! 			FLOW FIELD INITIATION           !
!===========================================!
U1 = 0.0D0
DO J = 2, N - 1
	U1(1, J) = 4 * ET(J) * (1 - ET(J))
	U1(M, J) = 4 * ET(J) * (1 - ET(J))
END DO

V1 = 0.0D0
P1 = 0.0D0
P1(:, 1) = 1.0D0

U2 = U1
V2 = V1
P2 = P1

!===========================================!
! 			    TIME MARCHING               !
!===========================================!
RES  = 10000
ITER = 1

OPEN(99, FILE = 'RES.PLT')

DO WHILE (RES .GT. EPS)
	! INTERIOR POINTS
	DO I = 2, M - 1
		DO J = 2, N - 1					   
			OMEGA = 1 + 2 * dT / RE * ((dXdXI(I) ** 2) / (dX ** 2) + &
			        (dYdET(J) ** 2) / (dY ** 2))
			
			U3(I, J) = 1.0D0 / OMEGA * (U1(I, J) - dT * (dXdXI(I) * &
					   (U2(I + 1, J) ** 2 - U2(I - 1, J) ** 2) / dX &
					   + dYdET(J) * (U2(I, J + 1) * V2(I, J + 1) -  &
					   U2(I, J - 1) * V2(I, J - 1)) / dY + dXdXI(I) &
					   * (P2(I + 1, J) - P2(I - 1, J)) / dX) + 2 *  &
					   dT / RE * ((dXdXI(I) ** 2) * (U2(I + 1, J)   &
					   - U1(I, J) + U2(I - 1, J)) / (dX ** 2) +     &
					   (dYdET(J) ** 2) * (U2(I, J + 1) - U1(I, J) + &
					   U2(I, J - 1)) / (dY ** 2)))
					   
			V3(I, J) = 1.0D0 / OMEGA * (V1(I, J) - dT * (dXdXI(I) * &
					   (U2(I + 1, J) * V2(I + 1, J) - U2(I - 1, J)  &
					   * V2(I - 1, J)) / dX + dYdET(J) * (V2(I, J + &
					   1) ** 2 -  V2(I, J - 1) ** 2) / dY + dYdET(J)&
					   * (P2(I, J + 1) - P2(I, J - 1)) / dY) + 2 *  &
					   dT / RE * ((dXdXI(I) ** 2) * (V2(I + 1, J)   &
					   - V1(I, J) + V2(I - 1, J)) / (dX ** 2) +     &
					   (dYdET(J) ** 2) * (V2(I, J + 1) - V1(I, J) + &
					   V2(I, J - 1)) / (dY ** 2)))
            
            P3(I, J) = P1(I, J) - dT / DELT * (dXdXI(I) * (U2(I + 1, J) &
                       - U2(I - 1, J)) / dX +  dYdET(J) * (V2(I, J + 1) &
                       - V2(I, J - 1)) / dY)
		END DO
	END DO

!===========================================!
! 		 SPECIFY BOUNDARY CONDITIONS        !
!===========================================!
	! LOWER BOUNDARY
	J = 1
    DO I = 1, M
        U3(I, J) = 0.0D0
		V3(I, J) = 0.0D0
    END DO
    
	DO I = 2, M - 1
		P3(I, J) = P1(I, J) - 2 * dT / DELT * (dXdXI(I) * (U2(I + 1, J) - &
		           U2(I - 1, J)) / (2 * dX) + dYdET(J) * (V2(I, J + 1) -  &
				   V2(I, J)) / dY)
	END DO
	
	! UPPER BOUNDARY
	J = N
    DO I = 1, M
		U3(I, J) = 0.0D0
		V3(I, J) = 0.0D0
    END DO
    
	DO I = 2, M - 1
		P3(I, J) = P1(I, J) - 2 * dT / DELT * (dXdXI(I) * (U2(I + 1, J) - &
		           U2(I - 1, J)) / (2 * dX) + dYdET(J) * (V2(I, J) - V2   &
				   (I, J - 1)) / dY)
	END DO

	! INLET BOUNDARY
	I = 1
	DO J = 2, N - 1
		U3(I, J) = 4 * ET(J) * (1 - ET(J))
		V3(I, J) = 0.0D0
		P3(I, J) = P1(I, J) - 2 * dT / DELT * (dXdXI(I) * (U2(I + 1, J) - &
		           U2(I, J)) / dX + dYdET(J) * (V2(I, J + 1) - V2(I, J -  &
				   1)) / (2 * dY))
	END DO	
	
	! OUTLET NOUNDARY
	I = M
	DO J = 2, N - 1
		U3(I, J) = 4 * ET(J) * (1 - ET(J))
		V3(I, J) = 0.0D0
		P3(I, J) = P1(I, J) - 2 * dT / DELT * (dXdXI(I) * (U2(I, J) - &
		           U2(I - 1, J)) / dX + dYdET(J) * (V2(I, J + 1) - V2 &
				   (I, J - 1)) / (2 * dY))
	END DO

!===========================================!
!            RESIDUAL CALCULATION           !
!===========================================!
	RES = -100000.0D0
	DO I = 1, M
		DO J = 1, N
			TEMP = ABS(U3(I, J) - U2(I, J))
			IF (RES .LT. TEMP) THEN
				RES = TEMP
			END IF
		END DO
	END DO
	WRITE( *, *) 'THE', ITER, 'TH, RESIDUAL ERROR:', RES
	WRITE(99, *) ITER, RES
	
	ITER = ITER + 1
	
	U1 = U2
	V1 = V2
	P1 = P2
	
	U2 = U3
	V2 = V3
	P2 = P3	
    
    IF (ITER .GT. MAXITER) THEN
        EXIT
    END IF
    
END DO
CLOSE(99)

READ(*, *)

!===========================================!
! 		       OUTPUT RESULTS               !
!===========================================!
OPEN(100, FILE = 'FIELD.PLT')
WRITE(100,*) 'title = parameter'
WRITE(100,*) 'variables = X, Y, P, U, V'
WRITE(100,*) 'zone f=point I=', N, ',J=', M
DO I = 1, M
	DO J = 1, N
		WRITE(100,*) XI(I), ET(J), P3(I, J), U3(I, J), V3(I, J)
	END DO
ENDDO
CLOSE(100)


OPEN (101, FILE = 'X.DAT')
DO J = 1, N
    WRITE(101, '(1000F20.12)') XI
ENDDO
CLOSE(101)

OPEN (102, FILE = 'Y.DAT')
DO I = 1, M
    WRITE(102, '(1000F20.12)') ET
ENDDO
CLOSE(102)

OPEN (103, FILE = 'U.DAT')
DO J = 1, N
    WRITE(103, '(1000F20.12)') U3(:, J)
ENDDO
CLOSE(103)

OPEN (104, FILE = 'V.DAT')
DO J = 1, N
    WRITE(104, '(1000F20.12)') V3(:, J)
ENDDO
CLOSE(104)

OPEN (105, FILE = 'P.DAT')
DO J = 1, N
    WRITE(105, '(1000F20.12)') P3(:, J)
ENDDO
CLOSE(105)

END PROGRAM
